Validation of a stereological method for estimating particle size and density from 2D projections with high accuracy

Stereological methods for estimating the 3D particle size and density from 2D projections are essential to many research fields. These methods are, however, prone to errors arising from undetected particle profiles due to sectioning and limited resolution, known as ‘lost caps’. A potential solution developed by Keiding, Jensen, and Ranek in 1972, which we refer to as the Keiding model, accounts for lost caps by quantifying the smallest detectable profile in terms of its limiting ‘cap angle’ (ϕ), a size-independent measure of a particle’s distance from the section surface. However, this simple solution has not been widely adopted nor tested. Rather, model-independent design-based stereological methods, which do not explicitly account for lost caps, have come to the fore. Here, we provide the first experimental validation of the Keiding model by comparing the size and density of particles estimated from 2D projections with direct measurement from 3D EM reconstructions of the same tissue. We applied the Keiding model to estimate the size and density of somata, nuclei and vesicles in the cerebellum of mice and rats, where high packing density can be problematic for design-based methods. Our analysis reveals a Gaussian distribution for ϕ rather than a single value. Nevertheless, curve fits of the Keiding model to the 2D diameter distribution accurately estimate the mean ϕ and 3D diameter distribution. While systematic testing using simulations revealed an upper limit to determining ϕ, our analysis shows that estimated ϕ can be used to determine the 3D particle density from the 2D density under a wide range of conditions, and this method is potentially more accurate than minimum-size-based lost-cap corrections and disector methods. Our results show the Keiding model provides an efficient means of accurately estimating the size and density of particles from 2D projections even under conditions of a high density.

2 Appendix S1. Derivation of Eq 3 and why dmin is not a good measure of lost caps. Abercrombie (1946) derived the following relation between the expected true particle count (Ntrue) and actual crude particle count (Nmeasured) of a section of thickness T, also known as the Abercrombie correction formula: Dividing both sides of this equation by the cross-sectional area over which the particles are counted (Areaxy) gives: Substituting in terms λ3D and λ2D used in this study gives: However, this relation is only correct if there are no lost caps. If there are lost caps, then μD must be scaled smaller. Using the Keiding fixed-ϕ model (1972), scaling of μD can be achieved via cosϕ (Fig 1) and λ3D can be estimated via the following (Cruz-Orive, 1983): 3 = 2 ( + cos ) Eq A1.2a which is Eq 3 of this study. Floderus (1944) derived a similar expression with respect to the section penetration depth of the smallest observable cap (hmin): and this expression was later recast by Konigsmark (1970)  1/2 Eq A1.5 Estimating particle size and density from 2D projections Rothman et al. 2023 3 For analysing images, Eq A1.4 has proved more useful than the original Floderus equation since dmin can be directly measured from the sample of measured 2D diameters. However, an underlying assumption of the Floderus/Konigsmark correction is that all particles are the same size. If there is a distribution of particle size, then each particle will have its own minimum observable diameter (δmin) and dmin will be the minimum of all δmin (Fig 6C1 and S8C1). Depending on the spread of the δmin distribution, the difference between dmin and the mean δmin can be large, in which case the dmin correction will create a large underestimation of λ3D. Moreover, there is a high probability that dmin is an outlier, i.e. an unusually small value; this could happen when conditions for identifying small caps are better than average, or when there are false positive measurements.
Our 3D analysis of MFT vesicles showed the Keiding model accurately describes the δmin-D relation of the vesicles (Fig 6C1 and S8C1) and therefore gives an accurate estimate of λ3D via Eq 3 (Fig 11). This is in contrast to the dmin correction which underestimated λ3D for the same datasets by 13-20% and the Abercrombie correction which underestimated λ3D by 22-24% ( Appendix S2. Estimation of the volume fraction of spherical particles from the area fraction of their 2D projections. The relation between the volume fraction (VF) of spherical particles and their observed area fraction (AF) in a 2D projection was derived by Weibel and Paumgartner (1978; their Eq 13 and 37) and is as follows: Here, m2 and m3 are dimensionless moments for a Gaussian distribution and μD = 2μR.
Because ϕ is a better descriptor of the lost-cap distribution than hmin or its equivalent dmin (Appendix S1; Fig 6C1 and S8C1), we express X as a function of ϕ rather than hmin by first defining X with respect to dmin via Eq A1.5: Eq A2.2b and substituting dϕ = μD·sinϕ for dmin: To test this relation, we can consider the ideal scenario of a planar section with no lost caps (ϕ = 0°), in which case g = 0, X = 0, Kv = 1 and VF = AF, which is expected (Weibel and Paumgartner, 1978). We can also consider the scenario where all particles have the same diameter (D), in which case σD = 0, m2 = m3 = 1 and:  G(d) of liver cell nuclei (blue circles) computed from Table 1 of Keiding et al. (1972) for patients 601, 2003 and 1037. Distribution x-axes radii (units of mm) were converted to diameters (units of μm) via the following conversion factor: Θ = 2×1000/2300 μm/mm, where 2300 corrects for magnification. The maximum likelihood estimation (MLE) fits of Keiding et al. (black lines) were computed by plugging the parameters from their Table 2 (f, β, ϕ, p1, p2) into a modified version of Eq 1 (NMKeidingChi3) that includes the weighted sum of 3 G(d) as described by Keiding et al. (their Eq 4.1) using a chi distribution for F(d) (K-Chi3; Eq 6) and converting β from units of square radii (mm 2 ) to square diameters (μm 2 ) by multiplying by Θ 2 . For comparison, MLE fits were recomputed using a Gaussian distribution for F(d) (red lines; K-Gauss3; NMKeidingGauss3; Eq 5) where μD and σD were computed from f and β. The overlap of these two curves (K-Chi3 vs. K-Gauss3) demonstrates the estimated chi distributions from the MLE fits are approximately Gaussian, which is expected since f is large for these fits (104, 208 and 70). As a last comparison, the same modified version of Eq 1 (NMKeidingGauss3) was curve fitted to the 3 G(d) using our LSE routine (blue dashed lines) resulting in nearly equivalent curves and parameters (Table  4).  Δϕ of Keiding-model fits to simulated G(d) versus the difference between estimated ϕ and estimated ϕcutoff, where estimated ϕcutoff was computed via Eq 9 using estimated μD and σD. These plots show that if estimated ϕ < estimated ϕcutoff (left side of graphs) there is a high probability |Δϕ| < 4° for G(d) computed from ~500 diameters (A1; red shading; data from Fig 4C), |Δϕ| < 3° for G(d) computed from ~2000 diameters (A2; data from Fig S4) and |Δϕ| < 5° for G(d) computed from ~300 diameters (not shown).

Fig S6. 3D analysis of MFT vesicles in ET z-stacks.
A. Example curve fit of an ellipse (Eq 10; red solid line) to the equivalent xy-radius (½darea) vs. z-stack image number (z#) relation of a MFT vesicle (black circles; vesicle #20) from z-stack ET10 where the z-axis was centered at z# = 0. Top graph shows residuals (nm) between data and fit. Red dashed lines denote measured ϕ for the vesicle's north and south poles computed as ϕ = sin -1 (δmin/D), where δmin is the minimum diameter measured at the given pole.

B.
Probability density of eccentricity factor (E; 0.1 bins) for curve fits of Eq 10 to the xyradius vs. z# relations, as in A, of vesicles in ET10 (red; n = 132) and ET11 (blue; n = 233). For ET10 and ET11, Sz was fixed to 0.63 and 0.53 nm, respectively, the values necessary to obtain an average E = 1.00, i.e. an average z-axis diameter equal to the average xy-axis diameter (isotropic conditions). These Sz are 1.7-and 1.4-fold larger than the acquisition Sz (0.38 nm) indicating the sections after imaging were 60% and 70% of their original thickness, a shrinkage that is consistent with previous findings (Luther et al. 1988  analysis; assuming true ϕ = 5°). Inset: cartoon of a particle sectioned within a z-stack.
Dashed red lines denote the particle's ϕ.

Fig S7. F(d) of cerebellar MFT vesicles and GC nuclei is well described by a Gaussian distribution.
Comparison of LSE curve fits of a Gaussian distribution (Eq 5; red line), chi distribution  Outlines were drawn around the outer contour of visually identified GC nuclei (yellow; blind detection) and darea computed from the outlines. A subset of nuclei (n = 30) were tracked through multiple z-planes in steps of 5 (0.2 μm) with one shown here (cyan; nonblind detection; overlaid outlines from images above and below). Pink dashed line illustrates a measurement of a nuclei's xy-orientation with respect to its long diameter (θlong = 144°). Only a subregion of the analysis is shown. Scale bar 6 μm.
B. Equivalent xy-radius (r = ½darea) vs. z# relation of the 30 nuclei tracked through multiple z-planes (lines), including the nucleus in A (cyan circles). The r-z# relations were curve fitted to an ellipse (Eq 10), with one fit shown here for the nucleus in A (black line; D = 6.47 ± 0.04 μm, E = 0.96 ± 0.01). For all fits, Sz was fixed to 40 nm, the value necessary to obtain an average E = 1.00, i.e. an average z-axis diameter approximately equal to the average xy-axis diameter (isotropic conditions). This Sz is similar to that reported by Nguyen et al. (45 nm). Z-axes were centered at z# = 0. Inset: probability density of E (black circles; 0.15 bins) curve fitted to a Gaussian function (gray line; Eq. 5). The narrow E distribution (σ = ±0.17) indicates spherical dimensions.
Measures of δmin and ϕ are not reported due to a large discretization error for this dataset (Fig S6D).
C. Probability density (per °) of θlong, illustrated in A, showing the GC nuclei have no systematic orientation in the xy-plane, consistent with a random orientation (n = 107; 20° bins). For a given nucleus, θlong was measured on the z-plane where the nucleus has a maximum darea, described in D.
Only a subset of the original TEM z-stack was analysed: coordinate pixels x = 124016-146544, y = 95976-108264, z = 70-495 in z-steps of 5; however, images for z = 135, 175 and 445 were non-existent in the online database and substituted with z = 136, 176 and 446. Images were down-sampled by a factor of 10, giving Sxy = 40 nm/pixel. C. Average ΔμD, ΔσD and Δϕ of curve fits of Eq 1 to the 100 G(d) of the best-match simulations in B1, computed with respect to 'measured' μD, σD and μϕ of the particles in the projection. Black dashed lines denote ΔμD, ΔσD and Δϕ of the experimental data.
Errors of the Gaussian-ϕ model match the experimental errors better than those of the fix-ϕ model.
To find the best-match μϕ for the given experimental G(d), the sum of squared differences (χ 2 ) was computed between the experimental G(d) and simulated G(d).
Cumulative distribution functions (CDFs) of χ 2 were computed from the 100 simulated G(d) at a given μϕ. From the CDFs, χ 2 values with 50% probability (χ 2 -50%) were computed as a function of μϕ, as well as χ 2 values with 16 and 84% probabilities representing ±σ. The μϕ with smallest χ 2 -50% was deemed the best match. Confidence intervals above and below the best-match μϕ denote the range of μϕ whose χ 2 are not significantly different to that of the best-match μϕ, computed via a KS test (p > 0.05).

Fig S11. The Keiding-model fixed-ϕ assumption introduces only small errors when estimating F(d) and ϕ from G(d) (2D projection simulations).
Average   measured from TEM images of the GC layer (Fig 2C1 and 2C2) where each row is for a different mouse and each column is for a different MFT. Each G(d) was curve fitted to Eq 1 (red solid lines) resulting in estimates of F(d) (red dotted lines and circles). For all fits, ϕ has a large error and is greater than estimated ϕcutoff (Eq 9) indicating ϕ is indeterminable and G(d) ≈ F(d) (Fig 4C). Gaussian fits to the same G(d) (Eq 5; black dashed lines and circles) overlap the Keiding-model fits and estimated F(d).